###
###
### Clark, Zucker, and Urpelainen; Renewable and Sustainable Energy Reviews
### The Future of Coal-Fired Power Generation in Southeast Asia
###
### Replication file: generating datasets for descriptive statistics, plots
### Created 10 Dec 2019
###
###

############# #
#### SETUP ####
############# #

rm(list=ls())

pacman::p_load(dplyr,
               gridExtra,
               xtable,
               ggplot2,
               ggthemes,
               grid,
               tidyverse)



###################################### #
#### **** STANDARD CALCULATION **** ####
###################################### #

######################### #
#### BASELINE SCENARIO ####
######################### #

# > Data loading ####

secoal_orig <- read_csv("SEAsia_CoalPlants.csv")
secoal_orig <- secoal_orig[secoal_orig$region == "SE Asia" & secoal_orig$status != "retired",]
secoal <- secoal_orig[,c("capacity","status","country","commission_year")]

secoal_cyear <- read.csv("coal_cyear_seasia.csv")

status_percentages <- readRDS("status_percentages.rds")

secoal$commission_year <- as.numeric(as.character(secoal$commission_year))

capfact <- data.frame(country = c("Vietnam","Indonesia","Malaysia","Philippines","Mean"),
                      low = c(53,51,51,55,NA), high = c(70,74,74,78,NA),
                      mid = c(61.5,62.5,58,66.5,NA))
capfact$low[5] <- mean(capfact$low[1:4])
capfact$high[5] <- mean(capfact$high[1:4])
capfact$mid[5] <- mean(capfact$mid[1:4])
capfact[,2:4] <- capfact[,2:4]/100

rownames(capfact) <- capfact$country

# > Setting parameters ####

# operational lifespan for individual units
lifespan <- 35

# year span for projections
projspan <- 2017:2057

# 0: region-wide attrition rate, 1: country-specific
attrition_specificity <- 1

# 1: middle attrition rate; 1<: high attrition rate; 1>: low attrition rate
attrition_scenario <- 1

# "low": low capacity factor; "mid": mid capacity factor; "high": high capacity factor
cap_scenario <- "mid"

# est. commissioning year for operating units w/o a listed commissioning year
operating_year <- round(mean(secoal$commission_year, na.rm = T), 0)

# pipeline times
announced_time <- 10
prepermit_time <- 7
permitted_time <- 4
construction_time <- 1

# > Calculating projections ####

projspan <- paste0("Y",projspan)
lastyear <- as.numeric(substr(last(projspan), 2, 5))
firstyear <- as.numeric(substr(first(projspan), 2, 5))

attr_All <- (status_percentages$All[3] + status_percentages$All[7])/
  (status_percentages$All[3] + status_percentages$All[7] + status_percentages$All[6])*attrition_scenario
attr_Cambodia <- (status_percentages$Cambodia[3] + status_percentages$Cambodia[7])/
  (status_percentages$Cambodia[3] + status_percentages$Cambodia[7] + status_percentages$Cambodia[6])*attrition_scenario
attr_Indonesia <- (status_percentages$Indonesia[3] + status_percentages$Indonesia[7])/
  (status_percentages$Indonesia[3] + status_percentages$Indonesia[7] + status_percentages$Indonesia[6])*attrition_scenario
attr_Laos <- (status_percentages$Laos[3] + status_percentages$Laos[7])/
  (status_percentages$Laos[3] + status_percentages$Laos[7] + status_percentages$Laos[6])*attrition_scenario
attr_Malaysia <- (status_percentages$Malaysia[3] + status_percentages$Malaysia[7])/
  (status_percentages$Malaysia[3] + status_percentages$Malaysia[7] + status_percentages$Malaysia[6])*attrition_scenario
attr_Myanmar <- (status_percentages$Myanmar[3] + status_percentages$Myanmar[7])/
  (status_percentages$Myanmar[3] + status_percentages$Myanmar[7] + status_percentages$Myanmar[6])*attrition_scenario
attr_Philippines <- (status_percentages$Philippines[3] + status_percentages$Philippines[7])/
  (status_percentages$Philippines[3] + status_percentages$Philippines[7] + status_percentages$Philippines[6])*attrition_scenario
attr_Thailand <- (status_percentages$Thailand[3] + status_percentages$Thailand[7])/
  (status_percentages$Thailand[3] + status_percentages$Thailand[7] + status_percentages$Thailand[6])*attrition_scenario
attr_Vietnam <- (status_percentages$Vietnam[3] + status_percentages$Vietnam[7])/
  (status_percentages$Vietnam[3] + status_percentages$Vietnam[7] + status_percentages$Vietnam[6])*attrition_scenario

attr_All <- ifelse(attr_All > 1, 1, attr_All)
attr_Cambodia <- ifelse(attr_Cambodia > 1, 1, attr_Cambodia)
attr_Indonesia <- ifelse(attr_Indonesia > 1, 1, attr_Indonesia)
attr_Laos <- ifelse(attr_Laos > 1, 1, attr_Laos)
attr_Malaysia <- ifelse(attr_Malaysia > 1, 1, attr_Malaysia)
attr_Myanmar <- ifelse(attr_Myanmar > 1, 1, attr_Myanmar)
attr_Philippines <- ifelse(attr_Philippines > 1, 1, attr_Philippines)
attr_Thailand <- ifelse(attr_Thailand > 1, 1, attr_Thailand)
attr_Vietnam <- ifelse(attr_Vietnam > 1, 1, attr_Vietnam)

if(attrition_specificity == 0){
  secoal$capacity[secoal$status != "operating" & !is.na(secoal$capacity)] <- secoal$capacity[
    secoal$status != "operating" & !is.na(secoal$capacity)]*(1-attr_All)
} else{
  for (i in 1:nrow(secoal)){
    secoal$capacity[i] <- ifelse(secoal$country[i] == "Cambodia" &
                                   secoal$status[i] != "operating",
                                 secoal$capacity[i]*(1-attr_Cambodia),
                                 ifelse(secoal$country[i] == "Indonesia" &
                                          secoal$status[i] != "operating",
                                        secoal$capacity[i]*(1-attr_Indonesia),
                                        ifelse(secoal$country[i] == "Laos" &
                                                 secoal$status[i] != "operating",
                                               secoal$capacity[i]*(1-attr_Laos),
                                               ifelse(secoal$country[i] == "Malaysia" &
                                                        secoal$status[i] != "operating",
                                                      secoal$capacity[i]*(1-attr_Malaysia),
                                                      ifelse(secoal$country[i] == "Myanmar" &
                                                               secoal$status[i] != "operating",
                                                             secoal$capacity[i]*(1-attr_Myanmar),
                                                             ifelse(secoal$country[i] == "Philippines" &
                                                                      secoal$status[i] != "operating",
                                                                    secoal$capacity[i]*(1-attr_Philippines),
                                                                    ifelse(secoal$country[i] == "Thailand" &
                                                                             secoal$status[i] != "operating",
                                                                           secoal$capacity[i]*(1-attr_Thailand),
                                                                           ifelse(secoal$country[i] == "Vietnam" &
                                                                                    secoal$status[i] != "operating",
                                                                                  secoal$capacity[i]*(1-attr_Vietnam),
                                                                                  secoal$capacity[i]))))))))
  }
}

secoal$commission_year[secoal$status == "announced"] <- 2017 + announced_time
secoal$commission_year[secoal$status == "pre-permit"] <- 2017 + prepermit_time
secoal$commission_year[secoal$status == "permitted"] <- 2017 + permitted_time
secoal$commission_year[secoal$status == "construction"] <- 2017 + construction_time
secoal$commission_year[secoal$status == "operating" & is.na(secoal$commission_year)] <- operating_year

secoal <- secoal[!is.na(secoal$commission_year),]

secoal$retirement_year <- as.numeric(as.character(secoal$commission_year)) + lifespan

secoal[projspan] <- NA

for (j in 6:ncol(secoal)){
  .year <- as.numeric(substr(colnames(secoal)[j],2,5))
  for (i in 1:nrow(secoal)){
    secoal[i,j] <- ifelse(.year >= secoal$commission_year[i] &
                            .year <= secoal$retirement_year[i],
                          1, 0)
  }
}

secoal[,c(6:ncol(secoal))] <- sapply(secoal[,c(6:ncol(secoal))],
                                     FUN = function(x){secoal$capacity * x})
secoal[secoal$country == "Cambodia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Cambodia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Indonesia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Indonesia",c(6:ncol(secoal))],
                                                                  FUN = function(x){capfact["Indonesia",cap_scenario] * x})
secoal[secoal$country == "Laos",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Laos",c(6:ncol(secoal))],
                                                             FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Malaysia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Malaysia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Malaysia",cap_scenario] * x})
secoal[secoal$country == "Myanmar",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Myanmar",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Philippines",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Philippines",c(6:ncol(secoal))],
                                                                    FUN = function(x){capfact["Philippines",cap_scenario] * x})
secoal[secoal$country == "Thailand",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Thailand",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Vietnam",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Vietnam",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Vietnam",cap_scenario] * x})

all_proj <- data.frame(year = as.numeric(substr(projspan,2,5)), capacity = colSums(secoal[,6:ncol(secoal)]))
cambodia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Cambodia",6:ncol(secoal)]))
indonesia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                             capacity = colSums(secoal[secoal$country == "Indonesia",6:ncol(secoal)]))
laos_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                        capacity = colSums(secoal[secoal$country == "Laos",6:ncol(secoal)]))
malaysia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Malaysia",6:ncol(secoal)]))
myanmar_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Myanmar",6:ncol(secoal)]))
philippines_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                               capacity = colSums(secoal[secoal$country == "Philippines",6:ncol(secoal)]))
thailand_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Thailand",6:ncol(secoal)]))
vietnam_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Vietnam",6:ncol(secoal)]))

projectiondf <- data.frame(Year = c(firstyear,
                                    round((lastyear - firstyear)/4, 0) + firstyear,2*round((lastyear - firstyear)/4, 0) + firstyear,3*round((lastyear - firstyear)/4, 0) + firstyear,lastyear),
                           All = c(all_proj$capacity[all_proj$year == firstyear],
                                   all_proj$capacity[all_proj$year == round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == lastyear]),
                           Indonesia = c(indonesia_proj$capacity[indonesia_proj$year == firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == lastyear]),
                           Malaysia = c(malaysia_proj$capacity[malaysia_proj$year == firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == lastyear]),
                           Philippines = c(philippines_proj$capacity[philippines_proj$year == firstyear],
                                           philippines_proj$capacity[philippines_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == lastyear]),
                           Vietnam = c(vietnam_proj$capacity[vietnam_proj$year == firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == lastyear]))

projectiondf[,c(2:6)] <- round(projectiondf[,c(2:6)], 0)

countrylist <- c("All","Cambodia","Indonesia","Laos","Malaysia","Myanmar","Philippines","Thailand","Vietnam")
attrlist <- c(attr_All,attr_Cambodia,attr_Indonesia,attr_Laos,attr_Malaysia,attr_Myanmar,
              attr_Philippines,attr_Thailand,attr_Vietnam)
attrtab <- data.frame(Country = countrylist, Attrition = attrlist)
attrtab$Attrition <- round(attrtab$Attrition, 2)


saveRDS(attrtab, "midattrtab.rds")
saveRDS(projectiondf, "midprojectiondf.rds")


# > Creating plots ####

all_plot <- ggplot(all_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "All")
cambodia_plot <- ggplot(cambodia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Cambodia")
indonesia_plot <- ggplot(indonesia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Indonesia")
laos_plot <- ggplot(laos_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Laos")
malaysia_plot <- ggplot(malaysia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Malaysia")
myanmar_plot <- ggplot(myanmar_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Myanmar")
philippines_plot <- ggplot(philippines_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Philippines")
thailand_plot <- ggplot(thailand_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Thailand")
vietnam_plot <- ggplot(vietnam_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Vietnam")

grid.arrange(all_plot,cambodia_plot,indonesia_plot,laos_plot,malaysia_plot,myanmar_plot,
             philippines_plot,thailand_plot,vietnam_plot,
             ncol = 3, nrow = 3)

g <- arrangeGrob(all_plot,cambodia_plot,indonesia_plot,laos_plot,malaysia_plot,myanmar_plot,
                 philippines_plot,thailand_plot,vietnam_plot, ncol=3, nrow=3)

ggsave(file="mid3x3.pdf", g)


all_proj$country <- "SE Asia"
cambodia_proj$country <- "Cambodia"
indonesia_proj$country <- "Indonesia"
laos_proj$country <- "Laos"
malaysia_proj$country <- "Malaysia"
myanmar_proj$country <- "Myanmar"
philippines_proj$country <- "Philippines"
thailand_proj$country <- "Thailand"
vietnam_proj$country <- "Vietnam"

bycountry <- rbind(all_proj, cambodia_proj, indonesia_proj, laos_proj, malaysia_proj, philippines_proj,
                   thailand_proj, vietnam_proj)

bycountry$label <- bycountry$country
bycountry$label[bycountry$year != max(bycountry$year)-round((max(bycountry$year)-min(bycountry$year))/1.3,0)] <- ""

( midline <- ggplot(bycountry, aes(x = year, y = capacity, group = country, lty = as.factor(country))) +
    geom_line() + xlab("Year") + ylab("Effective Capacity (MW)") +
    scale_linetype_discrete(name = "Country") + theme_minimal() + 
    theme(
      panel.grid.major.x = element_blank(),
      panel.grid.major.y = element_blank(),
      panel.grid.minor.x = element_blank(),
      panel.grid.minor.y = element_blank(),
      axis.line.x = element_line(color="black", size = 0.1),
      axis.line.y = element_line(color="black", size = 0.1)) +
    ylim(0, 80000)
)


saveRDS(midline, "midline.rds")


# > CO2 estimates ####

co2f <- 8.2245*1000 * 215.225 * 0.00000000397347

bycountry_co2 <- bycountry
bycountry_co2$capacity <- bycountry_co2$capacity * co2f
colnames(bycountry_co2)[2] <- "co2_yearly"

bycountry_co2$co2_cumulative <- bycountry_co2$co2_yearly

for (i in 2:nrow(bycountry_co2)){
  bycountry_co2$co2_cumulative[i] <- ifelse(bycountry_co2$country[i-1] == bycountry_co2$country[i],
                                            bycountry_co2$co2_cumulative[i-1] + bycountry_co2$co2_cumulative[i],
                                            bycountry_co2$co2_cumulative[i])
}


saveRDS(bycountry_co2, "midco2.rds")


#################### #
#### LOW SCENARIO ####
#################### #

rm(list=ls())

# > Data loading ####

secoal_orig <- read_csv("SEAsia_CoalPlants.csv")
secoal_orig <- secoal_orig[secoal_orig$region == "SE Asia" & secoal_orig$status != "retired",]
secoal <- secoal_orig[,c("capacity","status","country","commission_year")]

secoal_cyear <- read.csv("coal_cyear_seasia.csv")

status_percentages <- readRDS("status_percentages.rds")

secoal$commission_year <- as.numeric(as.character(secoal$commission_year))

capfact <- data.frame(country = c("Vietnam","Indonesia","Malaysia","Philippines","Mean"),
                      low = c(53,51,51,55,NA), high = c(70,74,74,78,NA),
                      mid = c(61.5,62.5,58,66.5,NA))
capfact$low[5] <- mean(capfact$low[1:4])
capfact$high[5] <- mean(capfact$high[1:4])
capfact$mid[5] <- mean(capfact$mid[1:4])
capfact[,2:4] <- capfact[,2:4]/100

rownames(capfact) <- capfact$country

# > Setting parameters ####

# operational lifespan for individual units
lifespan <- 30

# year span for projections
projspan <- 2017:2057

# 0: region-wide attrition rate, 1: country-specific
attrition_specificity <- 1

# 1: middle attrition rate; 1<: high attrition rate; 1>: low attrition rate
attrition_scenario <- 1.1

# "low": low capacity factor; "mid": mid capacity factor; "high": high capacity factor
cap_scenario <- "low"

# est. commissioning year for operating units w/o a listed commissioning year
operating_year <- round(mean(secoal$commission_year, na.rm = T), 0)

# pipeline times

announced_time <- 10
prepermit_time <- 7
permitted_time <- 4
construction_time <- 1

# > Calculating projections ####

## AUTOMATED CALCULATIONS:

projspan <- paste0("Y",projspan)
lastyear <- as.numeric(substr(last(projspan), 2, 5))
firstyear <- as.numeric(substr(first(projspan), 2, 5))

attr_All <- (status_percentages$All[3] + status_percentages$All[7])/
  (status_percentages$All[3] + status_percentages$All[7] + status_percentages$All[6])*attrition_scenario
attr_Cambodia <- (status_percentages$Cambodia[3] + status_percentages$Cambodia[7])/
  (status_percentages$Cambodia[3] + status_percentages$Cambodia[7] + status_percentages$Cambodia[6])*attrition_scenario
attr_Indonesia <- (status_percentages$Indonesia[3] + status_percentages$Indonesia[7])/
  (status_percentages$Indonesia[3] + status_percentages$Indonesia[7] + status_percentages$Indonesia[6])*attrition_scenario
attr_Laos <- (status_percentages$Laos[3] + status_percentages$Laos[7])/
  (status_percentages$Laos[3] + status_percentages$Laos[7] + status_percentages$Laos[6])*attrition_scenario
attr_Malaysia <- (status_percentages$Malaysia[3] + status_percentages$Malaysia[7])/
  (status_percentages$Malaysia[3] + status_percentages$Malaysia[7] + status_percentages$Malaysia[6])*attrition_scenario
attr_Myanmar <- (status_percentages$Myanmar[3] + status_percentages$Myanmar[7])/
  (status_percentages$Myanmar[3] + status_percentages$Myanmar[7] + status_percentages$Myanmar[6])*attrition_scenario
attr_Philippines <- (status_percentages$Philippines[3] + status_percentages$Philippines[7])/
  (status_percentages$Philippines[3] + status_percentages$Philippines[7] + status_percentages$Philippines[6])*attrition_scenario
attr_Thailand <- (status_percentages$Thailand[3] + status_percentages$Thailand[7])/
  (status_percentages$Thailand[3] + status_percentages$Thailand[7] + status_percentages$Thailand[6])*attrition_scenario
attr_Vietnam <- (status_percentages$Vietnam[3] + status_percentages$Vietnam[7])/
  (status_percentages$Vietnam[3] + status_percentages$Vietnam[7] + status_percentages$Vietnam[6])*attrition_scenario

attr_All <- ifelse(attr_All > 1, 1, attr_All)
attr_Cambodia <- ifelse(attr_Cambodia > 1, 1, attr_Cambodia)
attr_Indonesia <- ifelse(attr_Indonesia > 1, 1, attr_Indonesia)
attr_Laos <- ifelse(attr_Laos > 1, 1, attr_Laos)
attr_Malaysia <- ifelse(attr_Malaysia > 1, 1, attr_Malaysia)
attr_Myanmar <- ifelse(attr_Myanmar > 1, 1, attr_Myanmar)
attr_Philippines <- ifelse(attr_Philippines > 1, 1, attr_Philippines)
attr_Thailand <- ifelse(attr_Thailand > 1, 1, attr_Thailand)
attr_Vietnam <- ifelse(attr_Vietnam > 1, 1, attr_Vietnam)

if(attrition_specificity == 0){
  secoal$capacity[secoal$status != "operating" & !is.na(secoal$capacity)] <- secoal$capacity[
    secoal$status != "operating" & !is.na(secoal$capacity)]*(1-attr_All)
} else{
  for (i in 1:nrow(secoal)){
    secoal$capacity[i] <- ifelse(secoal$country[i] == "Cambodia" &
                                   secoal$status[i] != "operating",
                                 secoal$capacity[i]*(1-attr_Cambodia),
                                 ifelse(secoal$country[i] == "Indonesia" &
                                          secoal$status[i] != "operating",
                                        secoal$capacity[i]*(1-attr_Indonesia),
                                        ifelse(secoal$country[i] == "Laos" &
                                                 secoal$status[i] != "operating",
                                               secoal$capacity[i]*(1-attr_Laos),
                                               ifelse(secoal$country[i] == "Malaysia" &
                                                        secoal$status[i] != "operating",
                                                      secoal$capacity[i]*(1-attr_Malaysia),
                                                      ifelse(secoal$country[i] == "Myanmar" &
                                                               secoal$status[i] != "operating",
                                                             secoal$capacity[i]*(1-attr_Myanmar),
                                                             ifelse(secoal$country[i] == "Philippines" &
                                                                      secoal$status[i] != "operating",
                                                                    secoal$capacity[i]*(1-attr_Philippines),
                                                                    ifelse(secoal$country[i] == "Thailand" &
                                                                             secoal$status[i] != "operating",
                                                                           secoal$capacity[i]*(1-attr_Thailand),
                                                                           ifelse(secoal$country[i] == "Vietnam" &
                                                                                    secoal$status[i] != "operating",
                                                                                  secoal$capacity[i]*(1-attr_Vietnam),
                                                                                  secoal$capacity[i]))))))))
  }
}

secoal$commission_year[secoal$status == "announced"] <- 2017 + announced_time
secoal$commission_year[secoal$status == "pre-permit"] <- 2017 + prepermit_time
secoal$commission_year[secoal$status == "permitted"] <- 2017 + permitted_time
secoal$commission_year[secoal$status == "construction"] <- 2017 + construction_time
secoal$commission_year[secoal$status == "operating" & is.na(secoal$commission_year)] <- operating_year

secoal <- secoal[!is.na(secoal$commission_year),]

secoal$retirement_year <- as.numeric(as.character(secoal$commission_year)) + lifespan

secoal[projspan] <- NA

for (j in 6:ncol(secoal)){
  .year <- as.numeric(substr(colnames(secoal)[j],2,5))
  for (i in 1:nrow(secoal)){
    secoal[i,j] <- ifelse(.year >= secoal$commission_year[i] &
                            .year <= secoal$retirement_year[i],
                          1, 0)
  }
}

secoal[,c(6:ncol(secoal))] <- sapply(secoal[,c(6:ncol(secoal))],
                                     FUN = function(x){secoal$capacity * x})
secoal[secoal$country == "Cambodia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Cambodia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Indonesia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Indonesia",c(6:ncol(secoal))],
                                                                  FUN = function(x){capfact["Indonesia",cap_scenario] * x})
secoal[secoal$country == "Laos",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Laos",c(6:ncol(secoal))],
                                                             FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Malaysia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Malaysia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Malaysia",cap_scenario] * x})
secoal[secoal$country == "Myanmar",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Myanmar",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Philippines",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Philippines",c(6:ncol(secoal))],
                                                                    FUN = function(x){capfact["Philippines",cap_scenario] * x})
secoal[secoal$country == "Thailand",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Thailand",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Vietnam",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Vietnam",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Vietnam",cap_scenario] * x})

all_proj <- data.frame(year = as.numeric(substr(projspan,2,5)), capacity = colSums(secoal[,6:ncol(secoal)]))
cambodia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Cambodia",6:ncol(secoal)]))
indonesia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                             capacity = colSums(secoal[secoal$country == "Indonesia",6:ncol(secoal)]))
laos_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                        capacity = colSums(secoal[secoal$country == "Laos",6:ncol(secoal)]))
malaysia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Malaysia",6:ncol(secoal)]))
myanmar_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Myanmar",6:ncol(secoal)]))
philippines_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                               capacity = colSums(secoal[secoal$country == "Philippines",6:ncol(secoal)]))
thailand_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Thailand",6:ncol(secoal)]))
vietnam_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Vietnam",6:ncol(secoal)]))

projectiondf <- data.frame(Year = c(firstyear,
                                    round((lastyear - firstyear)/4, 0) + firstyear,2*round((lastyear - firstyear)/4, 0) + firstyear,3*round((lastyear - firstyear)/4, 0) + firstyear,lastyear),
                           All = c(all_proj$capacity[all_proj$year == firstyear],
                                   all_proj$capacity[all_proj$year == round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == lastyear]),
                           Indonesia = c(indonesia_proj$capacity[indonesia_proj$year == firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == lastyear]),
                           Malaysia = c(malaysia_proj$capacity[malaysia_proj$year == firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == lastyear]),
                           Philippines = c(philippines_proj$capacity[philippines_proj$year == firstyear],
                                           philippines_proj$capacity[philippines_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == lastyear]),
                           Vietnam = c(vietnam_proj$capacity[vietnam_proj$year == firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == lastyear]))

projectiondf[,c(2:6)] <- round(projectiondf[,c(2:6)], 0)

countrylist <- c("All","Cambodia","Indonesia","Laos","Malaysia","Myanmar","Philippines","Thailand","Vietnam")
attrlist <- c(attr_All,attr_Cambodia,attr_Indonesia,attr_Laos,attr_Malaysia,attr_Myanmar,
              attr_Philippines,attr_Thailand,attr_Vietnam)
attrtab <- data.frame(Country = countrylist, Attrition = attrlist)
attrtab$Attrition <- round(attrtab$Attrition, 2)


saveRDS(attrtab, "lowattrtab.rds")
saveRDS(projectiondf, "lowprojectiondf.rds")


# > Creating plots ####

all_plot <- ggplot(all_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "All")
cambodia_plot <- ggplot(cambodia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Cambodia")
indonesia_plot <- ggplot(indonesia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Indonesia")
laos_plot <- ggplot(laos_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Laos")
malaysia_plot <- ggplot(malaysia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Malaysia")
myanmar_plot <- ggplot(myanmar_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Myanmar")
philippines_plot <- ggplot(philippines_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Philippines")
thailand_plot <- ggplot(thailand_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Thailand")
vietnam_plot <- ggplot(vietnam_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Vietnam")

grid.arrange(all_plot,cambodia_plot,indonesia_plot,laos_plot,malaysia_plot,myanmar_plot,
             philippines_plot,thailand_plot,vietnam_plot,
             ncol = 3, nrow = 3)

g <- arrangeGrob(all_plot,cambodia_plot,indonesia_plot,laos_plot,malaysia_plot,myanmar_plot,
                 philippines_plot,thailand_plot,vietnam_plot, ncol=3, nrow=3)

ggsave(file="low3x3.pdf", g)


all_proj$country <- "SE Asia"
cambodia_proj$country <- "Cambodia"
indonesia_proj$country <- "Indonesia"
laos_proj$country <- "Laos"
malaysia_proj$country <- "Malaysia"
myanmar_proj$country <- "Myanmar"
philippines_proj$country <- "Philippines"
thailand_proj$country <- "Thailand"
vietnam_proj$country <- "Vietnam"

bycountry <- rbind(all_proj, cambodia_proj, indonesia_proj, laos_proj, malaysia_proj, philippines_proj,
                   thailand_proj, vietnam_proj)

bycountry$label <- bycountry$country
bycountry$label[bycountry$year != max(bycountry$year)-round((max(bycountry$year)-min(bycountry$year))/1.3,0)] <- ""

( lowline <- ggplot(bycountry, aes(x = year, y = capacity, group = country, lty = as.factor(country))) +
    geom_line() + xlab("Year") + ylab("Effective Capacity (MW)") +
    scale_linetype_discrete(name = "Country") + theme_minimal() + 
    theme(legend.position = "none",
          panel.grid.major.x = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          axis.line.x = element_line(color="black", size = 0.1),
          axis.line.y = element_line(color="black", size = 0.1)) +
    ylim(0, 80000)
)


saveRDS(lowline, "lowline.rds")


# > CO2 estimates ####

co2f <- 7.528*1000 * 205.3 * 0.00000000397347

bycountry_co2 <- bycountry
bycountry_co2$capacity <- bycountry_co2$capacity * co2f
colnames(bycountry_co2)[2] <- "co2_yearly"

bycountry_co2$co2_cumulative <- bycountry_co2$co2_yearly

for (i in 2:nrow(bycountry_co2)){
  bycountry_co2$co2_cumulative[i] <- ifelse(bycountry_co2$country[i-1] == bycountry_co2$country[i],
                                            bycountry_co2$co2_cumulative[i-1] + bycountry_co2$co2_cumulative[i],
                                            bycountry_co2$co2_cumulative[i])
}


saveRDS(bycountry_co2, "lowco2_heat_emfactor.rds")


##################### #
#### HIGH SCENARIO ####
##################### #

rm(list=ls())

# > Data loading ####

secoal_orig <- read_csv("SEAsia_CoalPlants.csv")
secoal_orig <- secoal_orig[secoal_orig$region == "SE Asia" & secoal_orig$status != "retired",]
secoal <- secoal_orig[,c("capacity","status","country","commission_year")]

secoal_cyear <- read.csv("coal_cyear_seasia.csv")

status_percentages <- readRDS("status_percentages.rds")

secoal$commission_year <- as.numeric(as.character(secoal$commission_year))

capfact <- data.frame(country = c("Vietnam","Indonesia","Malaysia","Philippines","Mean"),
                      low = c(53,51,51,55,NA), high = c(70,74,74,78,NA),
                      mid = c(61.5,62.5,58,66.5,NA))
capfact$low[5] <- mean(capfact$low[1:4])
capfact$high[5] <- mean(capfact$high[1:4])
capfact$mid[5] <- mean(capfact$mid[1:4])
capfact[,2:4] <- capfact[,2:4]/100

rownames(capfact) <- capfact$country

# > Setting parameters ####

# operational lifespan for individual units
lifespan <- 40

# year span for projections
projspan <- 2017:2057

# 0: region-wide attrition rate, 1: country-specific
attrition_specificity <- 1

# 1: middle attrition rate; 1<: high attrition rate; 1>: low attrition rate
attrition_scenario <- .9

# "low": low capacity factor; "mid": mid capacity factor; "high": high capacity factor
cap_scenario <- "high"

# est. commissioning year for operating units w/o a listed commissioning year
operating_year <- round(mean(secoal$commission_year, na.rm = T), 0)

announced_time <- 10
prepermit_time <- 7
permitted_time <- 4
construction_time <- 1

# > Calculating projections ####

## AUTOMATED CALCULATIONS:

projspan <- paste0("Y",projspan)
lastyear <- as.numeric(substr(last(projspan), 2, 5))
firstyear <- as.numeric(substr(first(projspan), 2, 5))

attr_All <- (status_percentages$All[3] + status_percentages$All[7])/
  (status_percentages$All[3] + status_percentages$All[7] + status_percentages$All[6])*attrition_scenario
attr_Cambodia <- (status_percentages$Cambodia[3] + status_percentages$Cambodia[7])/
  (status_percentages$Cambodia[3] + status_percentages$Cambodia[7] + status_percentages$Cambodia[6])*attrition_scenario
attr_Indonesia <- (status_percentages$Indonesia[3] + status_percentages$Indonesia[7])/
  (status_percentages$Indonesia[3] + status_percentages$Indonesia[7] + status_percentages$Indonesia[6])*attrition_scenario
attr_Laos <- (status_percentages$Laos[3] + status_percentages$Laos[7])/
  (status_percentages$Laos[3] + status_percentages$Laos[7] + status_percentages$Laos[6])*attrition_scenario
attr_Malaysia <- (status_percentages$Malaysia[3] + status_percentages$Malaysia[7])/
  (status_percentages$Malaysia[3] + status_percentages$Malaysia[7] + status_percentages$Malaysia[6])*attrition_scenario
attr_Myanmar <- (status_percentages$Myanmar[3] + status_percentages$Myanmar[7])/
  (status_percentages$Myanmar[3] + status_percentages$Myanmar[7] + status_percentages$Myanmar[6])*attrition_scenario
attr_Philippines <- (status_percentages$Philippines[3] + status_percentages$Philippines[7])/
  (status_percentages$Philippines[3] + status_percentages$Philippines[7] + status_percentages$Philippines[6])*attrition_scenario
attr_Thailand <- (status_percentages$Thailand[3] + status_percentages$Thailand[7])/
  (status_percentages$Thailand[3] + status_percentages$Thailand[7] + status_percentages$Thailand[6])*attrition_scenario
attr_Vietnam <- (status_percentages$Vietnam[3] + status_percentages$Vietnam[7])/
  (status_percentages$Vietnam[3] + status_percentages$Vietnam[7] + status_percentages$Vietnam[6])*attrition_scenario

attr_All <- ifelse(attr_All > 1, 1, attr_All)
attr_Cambodia <- ifelse(attr_Cambodia > 1, 1, attr_Cambodia)
attr_Indonesia <- ifelse(attr_Indonesia > 1, 1, attr_Indonesia)
attr_Laos <- ifelse(attr_Laos > 1, 1, attr_Laos)
attr_Malaysia <- ifelse(attr_Malaysia > 1, 1, attr_Malaysia)
attr_Myanmar <- ifelse(attr_Myanmar > 1, 1, attr_Myanmar)
attr_Philippines <- ifelse(attr_Philippines > 1, 1, attr_Philippines)
attr_Thailand <- ifelse(attr_Thailand > 1, 1, attr_Thailand)
attr_Vietnam <- ifelse(attr_Vietnam > 1, 1, attr_Vietnam)

if(attrition_specificity == 0){
  secoal$capacity[secoal$status != "operating" & !is.na(secoal$capacity)] <- secoal$capacity[
    secoal$status != "operating" & !is.na(secoal$capacity)]*(1-attr_All)
} else{
  for (i in 1:nrow(secoal)){
    secoal$capacity[i] <- ifelse(secoal$country[i] == "Cambodia" &
                                   secoal$status[i] != "operating",
                                 secoal$capacity[i]*(1-attr_Cambodia),
                                 ifelse(secoal$country[i] == "Indonesia" &
                                          secoal$status[i] != "operating",
                                        secoal$capacity[i]*(1-attr_Indonesia),
                                        ifelse(secoal$country[i] == "Laos" &
                                                 secoal$status[i] != "operating",
                                               secoal$capacity[i]*(1-attr_Laos),
                                               ifelse(secoal$country[i] == "Malaysia" &
                                                        secoal$status[i] != "operating",
                                                      secoal$capacity[i]*(1-attr_Malaysia),
                                                      ifelse(secoal$country[i] == "Myanmar" &
                                                               secoal$status[i] != "operating",
                                                             secoal$capacity[i]*(1-attr_Myanmar),
                                                             ifelse(secoal$country[i] == "Philippines" &
                                                                      secoal$status[i] != "operating",
                                                                    secoal$capacity[i]*(1-attr_Philippines),
                                                                    ifelse(secoal$country[i] == "Thailand" &
                                                                             secoal$status[i] != "operating",
                                                                           secoal$capacity[i]*(1-attr_Thailand),
                                                                           ifelse(secoal$country[i] == "Vietnam" &
                                                                                    secoal$status[i] != "operating",
                                                                                  secoal$capacity[i]*(1-attr_Vietnam),
                                                                                  secoal$capacity[i]))))))))
  }
}

secoal$commission_year[secoal$status == "announced"] <- 2017 + announced_time
secoal$commission_year[secoal$status == "pre-permit"] <- 2017 + prepermit_time
secoal$commission_year[secoal$status == "permitted"] <- 2017 + permitted_time
secoal$commission_year[secoal$status == "construction"] <- 2017 + construction_time
secoal$commission_year[secoal$status == "operating" & is.na(secoal$commission_year)] <- operating_year

secoal <- secoal[!is.na(secoal$commission_year),]

secoal$retirement_year <- as.numeric(as.character(secoal$commission_year)) + lifespan

secoal[projspan] <- NA

for (j in 6:ncol(secoal)){
  .year <- as.numeric(substr(colnames(secoal)[j],2,5))
  for (i in 1:nrow(secoal)){
    secoal[i,j] <- ifelse(.year >= secoal$commission_year[i] &
                            .year <= secoal$retirement_year[i],
                          1, 0)
  }
}

secoal[,c(6:ncol(secoal))] <- sapply(secoal[,c(6:ncol(secoal))],
                                     FUN = function(x){secoal$capacity * x})
secoal[secoal$country == "Cambodia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Cambodia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Indonesia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Indonesia",c(6:ncol(secoal))],
                                                                  FUN = function(x){capfact["Indonesia",cap_scenario] * x})
secoal[secoal$country == "Laos",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Laos",c(6:ncol(secoal))],
                                                             FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Malaysia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Malaysia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Malaysia",cap_scenario] * x})
secoal[secoal$country == "Myanmar",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Myanmar",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Philippines",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Philippines",c(6:ncol(secoal))],
                                                                    FUN = function(x){capfact["Philippines",cap_scenario] * x})
secoal[secoal$country == "Thailand",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Thailand",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Vietnam",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Vietnam",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Vietnam",cap_scenario] * x})

all_proj <- data.frame(year = as.numeric(substr(projspan,2,5)), capacity = colSums(secoal[,6:ncol(secoal)]))
cambodia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Cambodia",6:ncol(secoal)]))
indonesia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                             capacity = colSums(secoal[secoal$country == "Indonesia",6:ncol(secoal)]))
laos_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                        capacity = colSums(secoal[secoal$country == "Laos",6:ncol(secoal)]))
malaysia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Malaysia",6:ncol(secoal)]))
myanmar_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Myanmar",6:ncol(secoal)]))
philippines_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                               capacity = colSums(secoal[secoal$country == "Philippines",6:ncol(secoal)]))
thailand_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Thailand",6:ncol(secoal)]))
vietnam_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Vietnam",6:ncol(secoal)]))

projectiondf <- data.frame(Year = c(firstyear,
                                    round((lastyear - firstyear)/4, 0) + firstyear,2*round((lastyear - firstyear)/4, 0) + firstyear,3*round((lastyear - firstyear)/4, 0) + firstyear,lastyear),
                           All = c(all_proj$capacity[all_proj$year == firstyear],
                                   all_proj$capacity[all_proj$year == round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == lastyear]),
                           Indonesia = c(indonesia_proj$capacity[indonesia_proj$year == firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == lastyear]),
                           Malaysia = c(malaysia_proj$capacity[malaysia_proj$year == firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == lastyear]),
                           Philippines = c(philippines_proj$capacity[philippines_proj$year == firstyear],
                                           philippines_proj$capacity[philippines_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == lastyear]),
                           Vietnam = c(vietnam_proj$capacity[vietnam_proj$year == firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == lastyear]))

projectiondf[,c(2:6)] <- round(projectiondf[,c(2:6)], 0)

countrylist <- c("All","Cambodia","Indonesia","Laos","Malaysia","Myanmar","Philippines","Thailand","Vietnam")
attrlist <- c(attr_All,attr_Cambodia,attr_Indonesia,attr_Laos,attr_Malaysia,attr_Myanmar,
              attr_Philippines,attr_Thailand,attr_Vietnam)
attrtab <- data.frame(Country = countrylist, Attrition = attrlist)
attrtab$Attrition <- round(attrtab$Attrition, 2)


saveRDS(attrtab, "highattrtab.rds")
saveRDS(projectiondf, "highprojectiondf.rds")


# > Creating plots ####

## PLOTS:

all_plot <- ggplot(all_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "All")
cambodia_plot <- ggplot(cambodia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Cambodia")
indonesia_plot <- ggplot(indonesia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Indonesia")
laos_plot <- ggplot(laos_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Laos")
malaysia_plot <- ggplot(malaysia_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Malaysia")
myanmar_plot <- ggplot(myanmar_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Myanmar")
philippines_plot <- ggplot(philippines_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Philippines")
thailand_plot <- ggplot(thailand_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Thailand")
vietnam_plot <- ggplot(vietnam_proj, aes(year, capacity)) + geom_line() + labs(x = "Year", y = "Eff. Cap. (MW)", title = "Vietnam")

grid.arrange(all_plot,cambodia_plot,indonesia_plot,laos_plot,malaysia_plot,myanmar_plot,
             philippines_plot,thailand_plot,vietnam_plot,
             ncol = 3, nrow = 3)

g <- arrangeGrob(all_plot,cambodia_plot,indonesia_plot,laos_plot,malaysia_plot,myanmar_plot,
                 philippines_plot,thailand_plot,vietnam_plot, ncol=3, nrow=3) #generates g

ggsave(file="high3x3.pdf", g)


all_proj$country <- "SE Asia"
cambodia_proj$country <- "Cambodia"
indonesia_proj$country <- "Indonesia"
laos_proj$country <- "Laos"
malaysia_proj$country <- "Malaysia"
myanmar_proj$country <- "Myanmar"
philippines_proj$country <- "Philippines"
thailand_proj$country <- "Thailand"
vietnam_proj$country <- "Vietnam"

bycountry <- rbind(all_proj, cambodia_proj, indonesia_proj, laos_proj, malaysia_proj, philippines_proj,
                   thailand_proj, vietnam_proj)

bycountry$label <- bycountry$country
bycountry$label[bycountry$year != max(bycountry$year)-round((max(bycountry$year)-min(bycountry$year))/1.3,0)] <- ""

( highline <- ggplot(bycountry, aes(x = year, y = capacity, group = country, lty = as.factor(country))) +
    geom_line() + xlab("Year") + ylab("Effective Capacity (MW)") +
    scale_linetype_discrete(name = "Country") + theme_minimal() + 
    theme(legend.position = "none",
          panel.grid.major.x = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          axis.line.x = element_line(color="black", size = 0.1),
          axis.line.y = element_line(color="black", size = 0.1)) +
    ylim(0, 80000)
)


saveRDS(highline, "highline.rds")


# > CO2 estimates ####

co2f <- 8.921*1000 * 227.4 * 0.00000000397347

bycountry_co2 <- bycountry
bycountry_co2$capacity <- bycountry_co2$capacity * co2f
colnames(bycountry_co2)[2] <- "co2_yearly"

bycountry_co2$co2_cumulative <- bycountry_co2$co2_yearly

for (i in 2:nrow(bycountry_co2)){
  bycountry_co2$co2_cumulative[i] <- ifelse(bycountry_co2$country[i-1] == bycountry_co2$country[i],
                                            bycountry_co2$co2_cumulative[i-1] + bycountry_co2$co2_cumulative[i],
                                            bycountry_co2$co2_cumulative[i])
}


saveRDS(bycountry_co2, "highco2_heat_emfactor.rds")


######################################### #
#### **** ALTERNATIVE CALCULATION **** ####
######################################### #

######################### #
#### BASELINE SCENARIO ####
######################### #

rm(list=ls())

# > Data loading ####

secoal_orig <- read_csv("SEAsia_CoalPlants.csv")
secoal_orig <- secoal_orig[secoal_orig$region == "SE Asia" & secoal_orig$status != "retired",]
secoal <- secoal_orig[,c("capacity","status","country","commission_year")]

secoal_cyear <- read.csv("coal_cyear_seasia.csv")

secoal$commission_year <- as.numeric(as.character(secoal$commission_year))

capfact <- data.frame(country = c("Vietnam","Indonesia","Malaysia","Philippines","Mean"),
                      low = c(53,51,51,55,NA), high = c(70,74,74,78,NA),
                      mid = c(61.5,62.5,58,66.5,NA))
capfact$low[5] <- mean(capfact$low[1:4])
capfact$high[5] <- mean(capfact$high[1:4])
capfact$mid[5] <- mean(capfact$mid[1:4])
capfact[,2:4] <- capfact[,2:4]/100

rownames(capfact) <- capfact$country

# > Setting parameters ####

# operational lifespan for individual units
lifespan <- 35

# year span for projections
projspan <- 2017:2057

# 0: region-wide attrition rate, 1: country-specific
attrition_specificity <- 1

# 1: middle attrition rate; 1>: high attrition rate; 1<: low attrition rate
attrition_scenario <- 1

# "low": low capacity factor; "mid": mid capacity factor; "high": high capacity factor
cap_scenario <- "mid"

# est. commissioning year for operating units w/o a listed commissioning year
operating_year <- round(mean(secoal$commission_year, na.rm = T), 0)

# pipeline times

announced_time <- 10
prepermit_time <- 7
permitted_time <- 4
construction_time <- 1

# > Calculating projections ####

status_percentages <- data.frame(Status = c("operating","construction","cancelled","shelved"), All = NA)
unique_countries <- unique(secoal$country)[1:8]
status_percentages[unique_countries] <- NA

status_percentages$Cambodia[1] <- sum(secoal$capacity[secoal$country == "Cambodia" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Indonesia[1] <- sum(secoal$capacity[secoal$country == "Indonesia" &
                                                         secoal$status == "operating" &
                                                         secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                       na.rm = T)
status_percentages$Laos[1] <- sum(secoal$capacity[secoal$country == "Laos" &
                                                    secoal$status == "operating" &
                                                    secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                  na.rm = T)
status_percentages$Malaysia[1] <- sum(secoal$capacity[secoal$country == "Malaysia" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Myanmar[1] <- sum(secoal$capacity[secoal$country == "Myanmar" &
                                                       secoal$status == "operating" &
                                                       secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                     na.rm = T)
status_percentages$Philippines[1] <- sum(secoal$capacity[secoal$country == "Philippines" &
                                                           secoal$status == "operating" &
                                                           secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                         na.rm = T)
status_percentages$Thailand[1] <- sum(secoal$capacity[secoal$country == "Thailand" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Vietnam[1] <- sum(secoal$capacity[secoal$country == "Vietnam" &
                                                       secoal$status == "operating" &
                                                       secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                     na.rm = T)

status_percentages$Cambodia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Cambodia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Indonesia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Indonesia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Laos[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Laos" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Malaysia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Malaysia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Myanmar[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Myanmar" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Philippines[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Philippines" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Thailand[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Thailand" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Vietnam[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Vietnam" & secoal$status == x & !is.na(secoal$capacity)])
}
)

status_percentages$All <- rowSums(status_percentages[,c(3:10)])

## AUTOMATED CALCULATIONS:

projspan <- paste0("Y",projspan)
lastyear <- as.numeric(substr(last(projspan), 2, 5))
firstyear <- as.numeric(substr(first(projspan), 2, 5))

attr_All <- (status_percentages$All[1] + status_percentages$All[2])/
  (status_percentages$All[1] + status_percentages$All[2] + status_percentages$All[3] + status_percentages$All[4])*attrition_scenario
attr_Cambodia <- (status_percentages$Cambodia[1] + status_percentages$Cambodia[2])/
  (status_percentages$Cambodia[1] + status_percentages$Cambodia[2] + status_percentages$Cambodia[3] + status_percentages$Cambodia[4])*attrition_scenario
attr_Indonesia <- (status_percentages$Indonesia[1] + status_percentages$Indonesia[2])/
  (status_percentages$Indonesia[1] + status_percentages$Indonesia[2] + status_percentages$Indonesia[3] + status_percentages$Indonesia[4])*attrition_scenario
attr_Laos <- (status_percentages$Laos[1] + status_percentages$Laos[2])/
  (status_percentages$Laos[1] + status_percentages$Laos[2] + status_percentages$Laos[3] + status_percentages$Laos[4])*attrition_scenario
attr_Malaysia <- (status_percentages$Malaysia[1] + status_percentages$Malaysia[2])/
  (status_percentages$Malaysia[1] + status_percentages$Malaysia[2] + status_percentages$Malaysia[3] + status_percentages$Malaysia[4])*attrition_scenario
attr_Myanmar <- (status_percentages$Myanmar[1] + status_percentages$Myanmar[2])/
  (status_percentages$Myanmar[1] + status_percentages$Myanmar[2] + status_percentages$Myanmar[3] + status_percentages$Myanmar[4])*attrition_scenario
attr_Philippines <- (status_percentages$Philippines[1] + status_percentages$Philippines[2])/
  (status_percentages$Philippines[1] + status_percentages$Philippines[2] + status_percentages$Philippines[3] + status_percentages$Philippines[4])*attrition_scenario
attr_Thailand <- (status_percentages$Thailand[1] + status_percentages$Thailand[2])/
  (status_percentages$Thailand[1] + status_percentages$Thailand[2] + status_percentages$Thailand[3] + status_percentages$Thailand[4])*attrition_scenario
attr_Vietnam <- (status_percentages$Vietnam[1] + status_percentages$Vietnam[2])/
  (status_percentages$Vietnam[1] + status_percentages$Vietnam[2] + status_percentages$Vietnam[3] + status_percentages$Vietnam[4])*attrition_scenario

attr_All <- ifelse(attr_All > 1, 1, attr_All)
attr_Cambodia <- ifelse(attr_Cambodia > 1, 1, attr_Cambodia)
attr_Indonesia <- ifelse(attr_Indonesia > 1, 1, attr_Indonesia)
attr_Laos <- ifelse(attr_Laos > 1, 1, attr_Laos)
attr_Malaysia <- ifelse(attr_Malaysia > 1, 1, attr_Malaysia)
attr_Myanmar <- ifelse(attr_Myanmar > 1, 1, attr_Myanmar)
attr_Philippines <- ifelse(attr_Philippines > 1, 1, attr_Philippines)
attr_Thailand <- ifelse(attr_Thailand > 1, 1, attr_Thailand)
attr_Vietnam <- ifelse(attr_Vietnam > 1, 1, attr_Vietnam)

if(attrition_specificity == 0){
  secoal$capacity[secoal$status != "operating" & !is.na(secoal$capacity)] <- secoal$capacity[
    secoal$status != "operating" & !is.na(secoal$capacity)]*(1-attr_All)
} else{
  for (i in 1:nrow(secoal)){
    secoal$capacity[i] <- ifelse(secoal$country[i] == "Cambodia" &
                                   secoal$status[i] != "operating",
                                 secoal$capacity[i]*(1-attr_Cambodia),
                                 ifelse(secoal$country[i] == "Indonesia" &
                                          secoal$status[i] != "operating",
                                        secoal$capacity[i]*(1-attr_Indonesia),
                                        ifelse(secoal$country[i] == "Laos" &
                                                 secoal$status[i] != "operating",
                                               secoal$capacity[i]*(1-attr_Laos),
                                               ifelse(secoal$country[i] == "Malaysia" &
                                                        secoal$status[i] != "operating",
                                                      secoal$capacity[i]*(1-attr_Malaysia),
                                                      ifelse(secoal$country[i] == "Myanmar" &
                                                               secoal$status[i] != "operating",
                                                             secoal$capacity[i]*(1-attr_Myanmar),
                                                             ifelse(secoal$country[i] == "Philippines" &
                                                                      secoal$status[i] != "operating",
                                                                    secoal$capacity[i]*(1-attr_Philippines),
                                                                    ifelse(secoal$country[i] == "Thailand" &
                                                                             secoal$status[i] != "operating",
                                                                           secoal$capacity[i]*(1-attr_Thailand),
                                                                           ifelse(secoal$country[i] == "Vietnam" &
                                                                                    secoal$status[i] != "operating",
                                                                                  secoal$capacity[i]*(1-attr_Vietnam),
                                                                                  secoal$capacity[i]))))))))
  }
}

secoal$commission_year[secoal$status == "announced"] <- 2017 + announced_time
secoal$commission_year[secoal$status == "pre-permit"] <- 2017 + prepermit_time
secoal$commission_year[secoal$status == "permitted"] <- 2017 + permitted_time
secoal$commission_year[secoal$status == "construction"] <- 2017 + construction_time
secoal$commission_year[secoal$status == "operating" & is.na(secoal$commission_year)] <- operating_year

secoal <- secoal[!is.na(secoal$commission_year),]

secoal$retirement_year <- as.numeric(as.character(secoal$commission_year)) + lifespan

secoal[projspan] <- NA

for (j in 6:ncol(secoal)){
  .year <- as.numeric(substr(colnames(secoal)[j],2,5))
  for (i in 1:nrow(secoal)){
    secoal[i,j] <- ifelse(.year >= secoal$commission_year[i] &
                            .year <= secoal$retirement_year[i],
                          1, 0)
  }
}

secoal[,c(6:ncol(secoal))] <- sapply(secoal[,c(6:ncol(secoal))],
                                     FUN = function(x){secoal$capacity * x})
secoal[secoal$country == "Cambodia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Cambodia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Indonesia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Indonesia",c(6:ncol(secoal))],
                                                                  FUN = function(x){capfact["Indonesia",cap_scenario] * x})
secoal[secoal$country == "Laos",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Laos",c(6:ncol(secoal))],
                                                             FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Malaysia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Malaysia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Malaysia",cap_scenario] * x})
secoal[secoal$country == "Myanmar",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Myanmar",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Philippines",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Philippines",c(6:ncol(secoal))],
                                                                    FUN = function(x){capfact["Philippines",cap_scenario] * x})
secoal[secoal$country == "Thailand",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Thailand",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Vietnam",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Vietnam",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Vietnam",cap_scenario] * x})

all_proj <- data.frame(year = as.numeric(substr(projspan,2,5)), capacity = colSums(secoal[,6:ncol(secoal)]))
cambodia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Cambodia",6:ncol(secoal)]))
indonesia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                             capacity = colSums(secoal[secoal$country == "Indonesia",6:ncol(secoal)]))
laos_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                        capacity = colSums(secoal[secoal$country == "Laos",6:ncol(secoal)]))
malaysia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Malaysia",6:ncol(secoal)]))
myanmar_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Myanmar",6:ncol(secoal)]))
philippines_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                               capacity = colSums(secoal[secoal$country == "Philippines",6:ncol(secoal)]))
thailand_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Thailand",6:ncol(secoal)]))
vietnam_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Vietnam",6:ncol(secoal)]))

projectiondf <- data.frame(Year = c(firstyear,
                                    round((lastyear - firstyear)/4, 0) + firstyear,2*round((lastyear - firstyear)/4, 0) + firstyear,3*round((lastyear - firstyear)/4, 0) + firstyear,lastyear),
                           All = c(all_proj$capacity[all_proj$year == firstyear],
                                   all_proj$capacity[all_proj$year == round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == lastyear]),
                           Indonesia = c(indonesia_proj$capacity[indonesia_proj$year == firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == lastyear]),
                           Malaysia = c(malaysia_proj$capacity[malaysia_proj$year == firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == lastyear]),
                           Philippines = c(philippines_proj$capacity[philippines_proj$year == firstyear],
                                           philippines_proj$capacity[philippines_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == lastyear]),
                           Vietnam = c(vietnam_proj$capacity[vietnam_proj$year == firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == lastyear]))

projectiondf[,c(2:6)] <- round(projectiondf[,c(2:6)], 0)

countrylist <- c("All","Cambodia","Indonesia","Laos","Malaysia","Myanmar","Philippines","Thailand","Vietnam")
attrlist <- c(attr_All,attr_Cambodia,attr_Indonesia,attr_Laos,attr_Malaysia,attr_Myanmar,
              attr_Philippines,attr_Thailand,attr_Vietnam)
attrtab <- data.frame(Country = countrylist, Attrition = attrlist)
attrtab$Attrition <- round(attrtab$Attrition, 2)


saveRDS(attrtab, "midattrtab_Calc2.rds")
saveRDS(projectiondf, "midprojectiondf_Calc2.rds")


# > Creating plots ####

all_proj$country <- "SE Asia"
cambodia_proj$country <- "Cambodia"
indonesia_proj$country <- "Indonesia"
laos_proj$country <- "Laos"
malaysia_proj$country <- "Malaysia"
myanmar_proj$country <- "Myanmar"
philippines_proj$country <- "Philippines"
thailand_proj$country <- "Thailand"
vietnam_proj$country <- "Vietnam"

bycountry <- rbind(all_proj, cambodia_proj, indonesia_proj, laos_proj, malaysia_proj, philippines_proj,
                   thailand_proj, vietnam_proj)

bycountry$label <- bycountry$country
bycountry$label[bycountry$year != max(bycountry$year)-round((max(bycountry$year)-min(bycountry$year))/1.3,0)] <- ""

( midline <- ggplot(bycountry, aes(x = year, y = capacity, group = country, linetype = as.factor(country))) +
    geom_line() + xlab("Year") + ylab("Capacity") +
    scale_linetype_discrete(name = "Country") + theme_minimal() + 
    ylim(0, 100000)
)


saveRDS(midline, "midline_Calc2.rds")


# > CO2 estimates ####

#################### #
#### LOW SCENARIO ####
#################### #

rm(list=ls())

# > Data loading ####

secoal_orig <- read_csv("SEAsia_CoalPlants.csv")
secoal_orig <- secoal_orig[secoal_orig$region == "SE Asia" & secoal_orig$status != "retired",]
secoal <- secoal_orig[,c("capacity","status","country","commission_year")]

secoal_cyear <- read.csv("coal_cyear_seasia.csv")

secoal$commission_year <- as.numeric(as.character(secoal$commission_year))

capfact <- data.frame(country = c("Vietnam","Indonesia","Malaysia","Philippines","Mean"),
                      low = c(53,51,51,55,NA), high = c(70,74,74,78,NA),
                      mid = c(61.5,62.5,58,66.5,NA))
capfact$low[5] <- mean(capfact$low[1:4])
capfact$high[5] <- mean(capfact$high[1:4])
capfact$mid[5] <- mean(capfact$mid[1:4])
capfact[,2:4] <- capfact[,2:4]/100

rownames(capfact) <- capfact$country

# > Setting parameters ####

# operational lifespan for individual units
lifespan <- 30

# year span for projections
projspan <- 2017:2057

# 0: region-wide attrition rate, 1: country-specific
attrition_specificity <- 1

# 1: middle attrition rate; 1>: high attrition rate; 1<: low attrition rate
attrition_scenario <- 1.1

# "low": low capacity factor; "mid": mid capacity factor; "high": high capacity factor
cap_scenario <- "low"

# est. commissioning year for operating units w/o a listed commissioning year
operating_year <- round(mean(secoal$commission_year, na.rm = T), 0)

# pipeline times

announced_time <- 10
prepermit_time <- 7
permitted_time <- 4
construction_time <- 1

# > Calculating projections ####

status_percentages <- data.frame(Status = c("operating","construction","cancelled","shelved"), All = NA)
unique_countries <- unique(secoal$country)[1:8]
status_percentages[unique_countries] <- NA

status_percentages$Cambodia[1] <- sum(secoal$capacity[secoal$country == "Cambodia" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Indonesia[1] <- sum(secoal$capacity[secoal$country == "Indonesia" &
                                                         secoal$status == "operating" &
                                                         secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                       na.rm = T)
status_percentages$Laos[1] <- sum(secoal$capacity[secoal$country == "Laos" &
                                                    secoal$status == "operating" &
                                                    secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                  na.rm = T)
status_percentages$Malaysia[1] <- sum(secoal$capacity[secoal$country == "Malaysia" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Myanmar[1] <- sum(secoal$capacity[secoal$country == "Myanmar" &
                                                       secoal$status == "operating" &
                                                       secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                     na.rm = T)
status_percentages$Philippines[1] <- sum(secoal$capacity[secoal$country == "Philippines" &
                                                           secoal$status == "operating" &
                                                           secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                         na.rm = T)
status_percentages$Thailand[1] <- sum(secoal$capacity[secoal$country == "Thailand" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Vietnam[1] <- sum(secoal$capacity[secoal$country == "Vietnam" &
                                                       secoal$status == "operating" &
                                                       secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                     na.rm = T)

status_percentages$Cambodia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Cambodia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Indonesia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Indonesia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Laos[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Laos" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Malaysia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Malaysia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Myanmar[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Myanmar" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Philippines[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Philippines" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Thailand[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Thailand" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Vietnam[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Vietnam" & secoal$status == x & !is.na(secoal$capacity)])
}
)

status_percentages$All <- rowSums(status_percentages[,c(3:10)])

## AUTOMATED CALCULATIONS:

projspan <- paste0("Y",projspan)
lastyear <- as.numeric(substr(last(projspan), 2, 5))
firstyear <- as.numeric(substr(first(projspan), 2, 5))

attr_All <- (status_percentages$All[1] + status_percentages$All[2])/
  (status_percentages$All[1] + status_percentages$All[2] + status_percentages$All[3] + status_percentages$All[4])*attrition_scenario
attr_Cambodia <- (status_percentages$Cambodia[1] + status_percentages$Cambodia[2])/
  (status_percentages$Cambodia[1] + status_percentages$Cambodia[2] + status_percentages$Cambodia[3] + status_percentages$Cambodia[4])*attrition_scenario
attr_Indonesia <- (status_percentages$Indonesia[1] + status_percentages$Indonesia[2])/
  (status_percentages$Indonesia[1] + status_percentages$Indonesia[2] + status_percentages$Indonesia[3] + status_percentages$Indonesia[4])*attrition_scenario
attr_Laos <- (status_percentages$Laos[1] + status_percentages$Laos[2])/
  (status_percentages$Laos[1] + status_percentages$Laos[2] + status_percentages$Laos[3] + status_percentages$Laos[4])*attrition_scenario
attr_Malaysia <- (status_percentages$Malaysia[1] + status_percentages$Malaysia[2])/
  (status_percentages$Malaysia[1] + status_percentages$Malaysia[2] + status_percentages$Malaysia[3] + status_percentages$Malaysia[4])*attrition_scenario
attr_Myanmar <- (status_percentages$Myanmar[1] + status_percentages$Myanmar[2])/
  (status_percentages$Myanmar[1] + status_percentages$Myanmar[2] + status_percentages$Myanmar[3] + status_percentages$Myanmar[4])*attrition_scenario
attr_Philippines <- (status_percentages$Philippines[1] + status_percentages$Philippines[2])/
  (status_percentages$Philippines[1] + status_percentages$Philippines[2] + status_percentages$Philippines[3] + status_percentages$Philippines[4])*attrition_scenario
attr_Thailand <- (status_percentages$Thailand[1] + status_percentages$Thailand[2])/
  (status_percentages$Thailand[1] + status_percentages$Thailand[2] + status_percentages$Thailand[3] + status_percentages$Thailand[4])*attrition_scenario
attr_Vietnam <- (status_percentages$Vietnam[1] + status_percentages$Vietnam[2])/
  (status_percentages$Vietnam[1] + status_percentages$Vietnam[2] + status_percentages$Vietnam[3] + status_percentages$Vietnam[4])*attrition_scenario

attr_All <- ifelse(attr_All > 1, 1, attr_All)
attr_Cambodia <- ifelse(attr_Cambodia > 1, 1, attr_Cambodia)
attr_Indonesia <- ifelse(attr_Indonesia > 1, 1, attr_Indonesia)
attr_Laos <- ifelse(attr_Laos > 1, 1, attr_Laos)
attr_Malaysia <- ifelse(attr_Malaysia > 1, 1, attr_Malaysia)
attr_Myanmar <- ifelse(attr_Myanmar > 1, 1, attr_Myanmar)
attr_Philippines <- ifelse(attr_Philippines > 1, 1, attr_Philippines)
attr_Thailand <- ifelse(attr_Thailand > 1, 1, attr_Thailand)
attr_Vietnam <- ifelse(attr_Vietnam > 1, 1, attr_Vietnam)

if(attrition_specificity == 0){
  secoal$capacity[secoal$status != "operating" & !is.na(secoal$capacity)] <- secoal$capacity[
    secoal$status != "operating" & !is.na(secoal$capacity)]*(1-attr_All)
} else{
  for (i in 1:nrow(secoal)){
    secoal$capacity[i] <- ifelse(secoal$country[i] == "Cambodia" &
                                   secoal$status[i] != "operating",
                                 secoal$capacity[i]*(1-attr_Cambodia),
                                 ifelse(secoal$country[i] == "Indonesia" &
                                          secoal$status[i] != "operating",
                                        secoal$capacity[i]*(1-attr_Indonesia),
                                        ifelse(secoal$country[i] == "Laos" &
                                                 secoal$status[i] != "operating",
                                               secoal$capacity[i]*(1-attr_Laos),
                                               ifelse(secoal$country[i] == "Malaysia" &
                                                        secoal$status[i] != "operating",
                                                      secoal$capacity[i]*(1-attr_Malaysia),
                                                      ifelse(secoal$country[i] == "Myanmar" &
                                                               secoal$status[i] != "operating",
                                                             secoal$capacity[i]*(1-attr_Myanmar),
                                                             ifelse(secoal$country[i] == "Philippines" &
                                                                      secoal$status[i] != "operating",
                                                                    secoal$capacity[i]*(1-attr_Philippines),
                                                                    ifelse(secoal$country[i] == "Thailand" &
                                                                             secoal$status[i] != "operating",
                                                                           secoal$capacity[i]*(1-attr_Thailand),
                                                                           ifelse(secoal$country[i] == "Vietnam" &
                                                                                    secoal$status[i] != "operating",
                                                                                  secoal$capacity[i]*(1-attr_Vietnam),
                                                                                  secoal$capacity[i]))))))))
  }
}

secoal$commission_year[secoal$status == "announced"] <- 2017 + announced_time
secoal$commission_year[secoal$status == "pre-permit"] <- 2017 + prepermit_time
secoal$commission_year[secoal$status == "permitted"] <- 2017 + permitted_time
secoal$commission_year[secoal$status == "construction"] <- 2017 + construction_time
secoal$commission_year[secoal$status == "operating" & is.na(secoal$commission_year)] <- operating_year

secoal <- secoal[!is.na(secoal$commission_year),]

secoal$retirement_year <- as.numeric(as.character(secoal$commission_year)) + lifespan

secoal[projspan] <- NA

for (j in 6:ncol(secoal)){
  .year <- as.numeric(substr(colnames(secoal)[j],2,5))
  for (i in 1:nrow(secoal)){
    secoal[i,j] <- ifelse(.year >= secoal$commission_year[i] &
                            .year <= secoal$retirement_year[i],
                          1, 0)
  }
}

secoal[,c(6:ncol(secoal))] <- sapply(secoal[,c(6:ncol(secoal))],
                                     FUN = function(x){secoal$capacity * x})
secoal[secoal$country == "Cambodia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Cambodia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Indonesia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Indonesia",c(6:ncol(secoal))],
                                                                  FUN = function(x){capfact["Indonesia",cap_scenario] * x})
secoal[secoal$country == "Laos",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Laos",c(6:ncol(secoal))],
                                                             FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Malaysia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Malaysia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Malaysia",cap_scenario] * x})
secoal[secoal$country == "Myanmar",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Myanmar",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Philippines",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Philippines",c(6:ncol(secoal))],
                                                                    FUN = function(x){capfact["Philippines",cap_scenario] * x})
secoal[secoal$country == "Thailand",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Thailand",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Vietnam",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Vietnam",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Vietnam",cap_scenario] * x})

all_proj <- data.frame(year = as.numeric(substr(projspan,2,5)), capacity = colSums(secoal[,6:ncol(secoal)]))
cambodia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Cambodia",6:ncol(secoal)]))
indonesia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                             capacity = colSums(secoal[secoal$country == "Indonesia",6:ncol(secoal)]))
laos_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                        capacity = colSums(secoal[secoal$country == "Laos",6:ncol(secoal)]))
malaysia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Malaysia",6:ncol(secoal)]))
myanmar_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Myanmar",6:ncol(secoal)]))
philippines_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                               capacity = colSums(secoal[secoal$country == "Philippines",6:ncol(secoal)]))
thailand_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Thailand",6:ncol(secoal)]))
vietnam_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Vietnam",6:ncol(secoal)]))

projectiondf <- data.frame(Year = c(firstyear,
                                    round((lastyear - firstyear)/4, 0) + firstyear,2*round((lastyear - firstyear)/4, 0) + firstyear,3*round((lastyear - firstyear)/4, 0) + firstyear,lastyear),
                           All = c(all_proj$capacity[all_proj$year == firstyear],
                                   all_proj$capacity[all_proj$year == round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == lastyear]),
                           Indonesia = c(indonesia_proj$capacity[indonesia_proj$year == firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == lastyear]),
                           Malaysia = c(malaysia_proj$capacity[malaysia_proj$year == firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == lastyear]),
                           Philippines = c(philippines_proj$capacity[philippines_proj$year == firstyear],
                                           philippines_proj$capacity[philippines_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == lastyear]),
                           Vietnam = c(vietnam_proj$capacity[vietnam_proj$year == firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == lastyear]))

projectiondf[,c(2:6)] <- round(projectiondf[,c(2:6)], 0)

countrylist <- c("All","Cambodia","Indonesia","Laos","Malaysia","Myanmar","Philippines","Thailand","Vietnam")
attrlist <- c(attr_All,attr_Cambodia,attr_Indonesia,attr_Laos,attr_Malaysia,attr_Myanmar,
              attr_Philippines,attr_Thailand,attr_Vietnam)
attrtab <- data.frame(Country = countrylist, Attrition = attrlist)
attrtab$Attrition <- round(attrtab$Attrition, 2)


saveRDS(attrtab, "lowattrtab_Calc2.rds")
saveRDS(projectiondf, "lowprojectiondf_Calc2.rds")


# > Creating plots ####

all_proj$country <- "SE Asia"
cambodia_proj$country <- "Cambodia"
indonesia_proj$country <- "Indonesia"
laos_proj$country <- "Laos"
malaysia_proj$country <- "Malaysia"
myanmar_proj$country <- "Myanmar"
philippines_proj$country <- "Philippines"
thailand_proj$country <- "Thailand"
vietnam_proj$country <- "Vietnam"

bycountry <- rbind(all_proj, cambodia_proj, indonesia_proj, laos_proj, malaysia_proj, philippines_proj,
                   thailand_proj, vietnam_proj)

bycountry$label <- bycountry$country
bycountry$label[bycountry$year != max(bycountry$year)-round((max(bycountry$year)-min(bycountry$year))/1.3,0)] <- ""

( lowline <- ggplot(bycountry, aes(x = year, y = capacity, group = country, linetype = as.factor(country))) +
    geom_line() + xlab("Year") + ylab("Capacity") +
    scale_linetype_discrete(name = "Country") + theme_minimal() + 
    ylim(0, 100000)
)


saveRDS(lowline, "lowline_Calc2.rds")


# > CO2 estimates ####

##################### #
#### HIGH SCENARIO ####
##################### #

rm(list=ls())

# > Data loading ####

secoal_orig <- read_csv("SEAsia_CoalPlants.csv")
secoal_orig <- secoal_orig[secoal_orig$region == "SE Asia" & secoal_orig$status != "retired",]
secoal <- secoal_orig[,c("capacity","status","country","commission_year")]

secoal_cyear <- read.csv("coal_cyear_seasia.csv")

secoal$commission_year <- as.numeric(as.character(secoal$commission_year))

capfact <- data.frame(country = c("Vietnam","Indonesia","Malaysia","Philippines","Mean"),
                      low = c(53,51,51,55,NA), high = c(70,74,74,78,NA),
                      mid = c(61.5,62.5,58,66.5,NA))
capfact$low[5] <- mean(capfact$low[1:4])
capfact$high[5] <- mean(capfact$high[1:4])
capfact$mid[5] <- mean(capfact$mid[1:4])
capfact[,2:4] <- capfact[,2:4]/100

rownames(capfact) <- capfact$country

# > Setting parameters ####

# operational lifespan for individual units
lifespan <- 40

# year span for projections
projspan <- 2017:2057

# 0: region-wide attrition rate, 1: country-specific
attrition_specificity <- 1

# 1: middle attrition rate; 1>: high attrition rate; 1<: low attrition rate
attrition_scenario <- .9

# "low": low capacity factor; "mid": mid capacity factor; "high": high capacity factor
cap_scenario <- "high"

# est. commissioning year for operating units w/o a listed commissioning year
operating_year <- round(mean(secoal$commission_year, na.rm = T), 0)

# pipeline times

announced_time <- 10
prepermit_time <- 7
permitted_time <- 4
construction_time <- 1

# > Calculating projections ####

status_percentages <- data.frame(Status = c("operating","construction","cancelled","shelved"), All = NA)
unique_countries <- unique(secoal$country)[1:8]
status_percentages[unique_countries] <- NA

status_percentages$Cambodia[1] <- sum(secoal$capacity[secoal$country == "Cambodia" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Indonesia[1] <- sum(secoal$capacity[secoal$country == "Indonesia" &
                                                         secoal$status == "operating" &
                                                         secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                       na.rm = T)
status_percentages$Laos[1] <- sum(secoal$capacity[secoal$country == "Laos" &
                                                    secoal$status == "operating" &
                                                    secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                  na.rm = T)
status_percentages$Malaysia[1] <- sum(secoal$capacity[secoal$country == "Malaysia" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Myanmar[1] <- sum(secoal$capacity[secoal$country == "Myanmar" &
                                                       secoal$status == "operating" &
                                                       secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                     na.rm = T)
status_percentages$Philippines[1] <- sum(secoal$capacity[secoal$country == "Philippines" &
                                                           secoal$status == "operating" &
                                                           secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                         na.rm = T)
status_percentages$Thailand[1] <- sum(secoal$capacity[secoal$country == "Thailand" &
                                                        secoal$status == "operating" &
                                                        secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                      na.rm = T)
status_percentages$Vietnam[1] <- sum(secoal$capacity[secoal$country == "Vietnam" &
                                                       secoal$status == "operating" &
                                                       secoal$commission_year >= 2010 & secoal$commission_year <= 2017],
                                     na.rm = T)

status_percentages$Cambodia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Cambodia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Indonesia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Indonesia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Laos[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Laos" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Malaysia[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Malaysia" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Myanmar[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Myanmar" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Philippines[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Philippines" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Thailand[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Thailand" & secoal$status == x & !is.na(secoal$capacity)])
}
)
status_percentages$Vietnam[2:4] <- sapply(status_percentages$Status[2:4], FUN = function(x){
  sum(secoal$capacity[secoal$country == "Vietnam" & secoal$status == x & !is.na(secoal$capacity)])
}
)

status_percentages$All <- rowSums(status_percentages[,c(3:10)])

## AUTOMATED CALCULATIONS:

projspan <- paste0("Y",projspan)
lastyear <- as.numeric(substr(last(projspan), 2, 5))
firstyear <- as.numeric(substr(first(projspan), 2, 5))

attr_All <- (status_percentages$All[1] + status_percentages$All[2])/
  (status_percentages$All[1] + status_percentages$All[2] + status_percentages$All[3] + status_percentages$All[4])*attrition_scenario
attr_Cambodia <- (status_percentages$Cambodia[1] + status_percentages$Cambodia[2])/
  (status_percentages$Cambodia[1] + status_percentages$Cambodia[2] + status_percentages$Cambodia[3] + status_percentages$Cambodia[4])*attrition_scenario
attr_Indonesia <- (status_percentages$Indonesia[1] + status_percentages$Indonesia[2])/
  (status_percentages$Indonesia[1] + status_percentages$Indonesia[2] + status_percentages$Indonesia[3] + status_percentages$Indonesia[4])*attrition_scenario
attr_Laos <- (status_percentages$Laos[1] + status_percentages$Laos[2])/
  (status_percentages$Laos[1] + status_percentages$Laos[2] + status_percentages$Laos[3] + status_percentages$Laos[4])*attrition_scenario
attr_Malaysia <- (status_percentages$Malaysia[1] + status_percentages$Malaysia[2])/
  (status_percentages$Malaysia[1] + status_percentages$Malaysia[2] + status_percentages$Malaysia[3] + status_percentages$Malaysia[4])*attrition_scenario
attr_Myanmar <- (status_percentages$Myanmar[1] + status_percentages$Myanmar[2])/
  (status_percentages$Myanmar[1] + status_percentages$Myanmar[2] + status_percentages$Myanmar[3] + status_percentages$Myanmar[4])*attrition_scenario
attr_Philippines <- (status_percentages$Philippines[1] + status_percentages$Philippines[2])/
  (status_percentages$Philippines[1] + status_percentages$Philippines[2] + status_percentages$Philippines[3] + status_percentages$Philippines[4])*attrition_scenario
attr_Thailand <- (status_percentages$Thailand[1] + status_percentages$Thailand[2])/
  (status_percentages$Thailand[1] + status_percentages$Thailand[2] + status_percentages$Thailand[3] + status_percentages$Thailand[4])*attrition_scenario
attr_Vietnam <- (status_percentages$Vietnam[1] + status_percentages$Vietnam[2])/
  (status_percentages$Vietnam[1] + status_percentages$Vietnam[2] + status_percentages$Vietnam[3] + status_percentages$Vietnam[4])*attrition_scenario

attr_All <- ifelse(attr_All > 1, 1, attr_All)
attr_Cambodia <- ifelse(attr_Cambodia > 1, 1, attr_Cambodia)
attr_Indonesia <- ifelse(attr_Indonesia > 1, 1, attr_Indonesia)
attr_Laos <- ifelse(attr_Laos > 1, 1, attr_Laos)
attr_Malaysia <- ifelse(attr_Malaysia > 1, 1, attr_Malaysia)
attr_Myanmar <- ifelse(attr_Myanmar > 1, 1, attr_Myanmar)
attr_Philippines <- ifelse(attr_Philippines > 1, 1, attr_Philippines)
attr_Thailand <- ifelse(attr_Thailand > 1, 1, attr_Thailand)
attr_Vietnam <- ifelse(attr_Vietnam > 1, 1, attr_Vietnam)

if(attrition_specificity == 0){
  secoal$capacity[secoal$status != "operating" & !is.na(secoal$capacity)] <- secoal$capacity[
    secoal$status != "operating" & !is.na(secoal$capacity)]*(1-attr_All)
} else{
  for (i in 1:nrow(secoal)){
    secoal$capacity[i] <- ifelse(secoal$country[i] == "Cambodia" &
                                   secoal$status[i] != "operating",
                                 secoal$capacity[i]*(1-attr_Cambodia),
                                 ifelse(secoal$country[i] == "Indonesia" &
                                          secoal$status[i] != "operating",
                                        secoal$capacity[i]*(1-attr_Indonesia),
                                        ifelse(secoal$country[i] == "Laos" &
                                                 secoal$status[i] != "operating",
                                               secoal$capacity[i]*(1-attr_Laos),
                                               ifelse(secoal$country[i] == "Malaysia" &
                                                        secoal$status[i] != "operating",
                                                      secoal$capacity[i]*(1-attr_Malaysia),
                                                      ifelse(secoal$country[i] == "Myanmar" &
                                                               secoal$status[i] != "operating",
                                                             secoal$capacity[i]*(1-attr_Myanmar),
                                                             ifelse(secoal$country[i] == "Philippines" &
                                                                      secoal$status[i] != "operating",
                                                                    secoal$capacity[i]*(1-attr_Philippines),
                                                                    ifelse(secoal$country[i] == "Thailand" &
                                                                             secoal$status[i] != "operating",
                                                                           secoal$capacity[i]*(1-attr_Thailand),
                                                                           ifelse(secoal$country[i] == "Vietnam" &
                                                                                    secoal$status[i] != "operating",
                                                                                  secoal$capacity[i]*(1-attr_Vietnam),
                                                                                  secoal$capacity[i]))))))))
  }
}

secoal$commission_year[secoal$status == "announced"] <- 2017 + announced_time
secoal$commission_year[secoal$status == "pre-permit"] <- 2017 + prepermit_time
secoal$commission_year[secoal$status == "permitted"] <- 2017 + permitted_time
secoal$commission_year[secoal$status == "construction"] <- 2017 + construction_time
secoal$commission_year[secoal$status == "operating" & is.na(secoal$commission_year)] <- operating_year

secoal <- secoal[!is.na(secoal$commission_year),]

secoal$retirement_year <- as.numeric(as.character(secoal$commission_year)) + lifespan

secoal[projspan] <- NA

for (j in 6:ncol(secoal)){
  .year <- as.numeric(substr(colnames(secoal)[j],2,5))
  for (i in 1:nrow(secoal)){
    secoal[i,j] <- ifelse(.year >= secoal$commission_year[i] &
                            .year <= secoal$retirement_year[i],
                          1, 0)
  }
}

secoal[,c(6:ncol(secoal))] <- sapply(secoal[,c(6:ncol(secoal))],
                                     FUN = function(x){secoal$capacity * x})
secoal[secoal$country == "Cambodia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Cambodia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Indonesia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Indonesia",c(6:ncol(secoal))],
                                                                  FUN = function(x){capfact["Indonesia",cap_scenario] * x})
secoal[secoal$country == "Laos",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Laos",c(6:ncol(secoal))],
                                                             FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Malaysia",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Malaysia",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Malaysia",cap_scenario] * x})
secoal[secoal$country == "Myanmar",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Myanmar",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Philippines",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Philippines",c(6:ncol(secoal))],
                                                                    FUN = function(x){capfact["Philippines",cap_scenario] * x})
secoal[secoal$country == "Thailand",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Thailand",c(6:ncol(secoal))],
                                                                 FUN = function(x){capfact["Mean",cap_scenario] * x})
secoal[secoal$country == "Vietnam",c(6:ncol(secoal))] <- sapply(secoal[secoal$country == "Vietnam",c(6:ncol(secoal))],
                                                                FUN = function(x){capfact["Vietnam",cap_scenario] * x})

all_proj <- data.frame(year = as.numeric(substr(projspan,2,5)), capacity = colSums(secoal[,6:ncol(secoal)]))
cambodia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Cambodia",6:ncol(secoal)]))
indonesia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                             capacity = colSums(secoal[secoal$country == "Indonesia",6:ncol(secoal)]))
laos_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                        capacity = colSums(secoal[secoal$country == "Laos",6:ncol(secoal)]))
malaysia_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Malaysia",6:ncol(secoal)]))
myanmar_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Myanmar",6:ncol(secoal)]))
philippines_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                               capacity = colSums(secoal[secoal$country == "Philippines",6:ncol(secoal)]))
thailand_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                            capacity = colSums(secoal[secoal$country == "Thailand",6:ncol(secoal)]))
vietnam_proj <- data.frame(year = as.numeric(substr(projspan,2,5)),
                           capacity = colSums(secoal[secoal$country == "Vietnam",6:ncol(secoal)]))

projectiondf <- data.frame(Year = c(firstyear,
                                    round((lastyear - firstyear)/4, 0) + firstyear,2*round((lastyear - firstyear)/4, 0) + firstyear,3*round((lastyear - firstyear)/4, 0) + firstyear,lastyear),
                           All = c(all_proj$capacity[all_proj$year == firstyear],
                                   all_proj$capacity[all_proj$year == round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                       firstyear],
                                   all_proj$capacity[all_proj$year == lastyear]),
                           Indonesia = c(indonesia_proj$capacity[indonesia_proj$year == firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                   firstyear],
                                         indonesia_proj$capacity[indonesia_proj$year == lastyear]),
                           Malaysia = c(malaysia_proj$capacity[malaysia_proj$year == firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                 firstyear],
                                        malaysia_proj$capacity[malaysia_proj$year == lastyear]),
                           Philippines = c(philippines_proj$capacity[philippines_proj$year == firstyear],
                                           philippines_proj$capacity[philippines_proj$year == round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                                       firstyear],
                                           philippines_proj$capacity[philippines_proj$year == lastyear]),
                           Vietnam = c(vietnam_proj$capacity[vietnam_proj$year == firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 2*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == 3*round((lastyear - firstyear)/4, 0) +
                                                               firstyear],
                                       vietnam_proj$capacity[vietnam_proj$year == lastyear]))

projectiondf[,c(2:6)] <- round(projectiondf[,c(2:6)], 0)

countrylist <- c("All","Cambodia","Indonesia","Laos","Malaysia","Myanmar","Philippines","Thailand","Vietnam")
attrlist <- c(attr_All,attr_Cambodia,attr_Indonesia,attr_Laos,attr_Malaysia,attr_Myanmar,
              attr_Philippines,attr_Thailand,attr_Vietnam)
attrtab <- data.frame(Country = countrylist, Attrition = attrlist)
attrtab$Attrition <- round(attrtab$Attrition, 2)


saveRDS(attrtab, "highattrtab_Calc2.rds")
saveRDS(projectiondf, "highprojectiondf_Calc2.rds")


# > Creating plots ####

all_proj$country <- "SE Asia"
cambodia_proj$country <- "Cambodia"
indonesia_proj$country <- "Indonesia"
laos_proj$country <- "Laos"
malaysia_proj$country <- "Malaysia"
myanmar_proj$country <- "Myanmar"
philippines_proj$country <- "Philippines"
thailand_proj$country <- "Thailand"
vietnam_proj$country <- "Vietnam"

bycountry <- rbind(all_proj, cambodia_proj, indonesia_proj, laos_proj, malaysia_proj, philippines_proj,
                   thailand_proj, vietnam_proj)

bycountry$label <- bycountry$country
bycountry$label[bycountry$year != max(bycountry$year)-round((max(bycountry$year)-min(bycountry$year))/1.3,0)] <- ""

( highline <- ggplot(bycountry, aes(x = year, y = capacity, group = country, linetype = as.factor(country))) +
    geom_line() + xlab("Year") + ylab("Capacity") +
    scale_linetype_discrete(name = "Country") + theme_minimal() + 
    ylim(0, 100000)
)


saveRDS(highline, "highline_Calc2.rds")


# > CO2 estimates ####

################################### #
#### **** OVERALL DATA PREP **** ####
################################### #

rm(list=ls())

# > Loading datasets ####

highline <- readRDS("highline.rds")
midline <- readRDS("midline.rds")
lowline <- readRDS("lowline.rds")

highattrtab <- readRDS("highattrtab.rds")
midattrtab <- readRDS("midattrtab.rds")
lowattrtab <- readRDS("lowattrtab.rds")

highprojectiondf <- readRDS("highprojectiondf.rds")
midprojectiondf <- readRDS("midprojectiondf.rds")
lowprojectiondf <- readRDS("lowprojectiondf.rds")

# > Exporting line plots ####

pdf("lowline.pdf", width = 5, height = 5)
lowline
dev.off()

pdf("midline.pdf", width = 5, height = 5)
midline
dev.off()

pdf("highline.pdf", width = 6, height = 5)
highline
dev.off()

# > Combining tables ####

# Attrition table
attrtab <- left_join(lowattrtab, midattrtab, by = "Country")
attrtab <- left_join(attrtab, highattrtab, by = "Country")
colnames(attrtab)[2:4] <- c("Low","Middle","High")

sink("attrtab.tex")
print(xtable(attrtab), include.rownames = F, include.colnames = T, booktabs = T, floating = F)
sink()

# 2047 projection table
low47 <- lowprojectiondf[4,]
mid47 <- midprojectiondf[4,]
high47 <- highprojectiondf[4,]

tab47 <- rbind(low47, mid47, high47)
tab47$Year[1:3] <- c("Low","Mid","High")

sink("tab47.tex")
print(xtable(tab47, digits = c(0,0,0,0,0,0,0)), include.rownames = F, include.colnames = T, booktabs = T, floating = F)
sink()

# > Loading datasets (alternative calculation) ####

highline_Calc2 <- readRDS("highline_Calc2.rds")
midline_Calc2 <- readRDS("midline_Calc2.rds")
lowline_Calc2 <- readRDS("lowline_Calc2.rds")

highattrtab_Calc2 <- readRDS("highattrtab_Calc2.rds")
midattrtab_Calc2 <- readRDS("midattrtab_Calc2.rds")
lowattrtab_Calc2 <- readRDS("lowattrtab_Calc2.rds")

highprojectiondf_Calc2 <- readRDS("highprojectiondf_Calc2.rds")
midprojectiondf_Calc2 <- readRDS("midprojectiondf_Calc2.rds")
lowprojectiondf_Calc2 <- readRDS("lowprojectiondf_Calc2.rds")

# > Exporting line plots (alternative calculation) ####

pdf("lowline_Calc2.pdf", width = 5, height = 5)
lowline_Calc2
dev.off()

pdf("midline_Calc2.pdf", width = 5, height = 5)
midline_Calc2
dev.off()

pdf("highline_Calc2.pdf", width = 6, height = 5)
highline_Calc2
dev.off()

# > Combining tables (alternative calculation) ####

# Attrition table
attrtab_Calc2 <- left_join(lowattrtab_Calc2, midattrtab_Calc2, by = "Country")
attrtab_Calc2 <- left_join(attrtab_Calc2, highattrtab_Calc2, by = "Country")
colnames(attrtab_Calc2)[2:4] <- c("Low","Middle","High")

sink("attrtab_Calc2.tex")
print(xtable(attrtab_Calc2), include.rownames = F, include.colnames = T, booktabs = T, floating = F)
sink()

# 2047 projection table
low47_Calc2 <- lowprojectiondf_Calc2[4,]
mid47_Calc2 <- midprojectiondf_Calc2[4,]
high47_Calc2 <- highprojectiondf_Calc2[4,]

tab47_Calc2 <- rbind(low47_Calc2, mid47_Calc2, high47_Calc2)
tab47_Calc2$Year[1:3] <- c("Low","Mid","High")

sink("tab47_Calc2.tex")
print(xtable(tab47_Calc2, digits = c(0,0,0,0,0,0,0)), include.rownames = F, include.colnames = T, booktabs = T, floating = F)
sink()

# > CO2 calcuations ####

highco2 <- readRDS("highco2.rds")
midco2 <- readRDS("midco2.rds")
lowco2 <- readRDS("lowco2.rds")

highco2_2 <- readRDS("highco2_heat_emfactor.rds")
lowco2_2 <- readRDS("lowco2_heat_emfactor.rds")

highco2$co2_cumulative <- highco2$co2_cumulative/1000
midco2$co2_cumulative <- midco2$co2_cumulative/1000
lowco2$co2_cumulative <- lowco2$co2_cumulative/1000

highco2$co2_yearly <- highco2$co2_yearly/1000
midco2$co2_yearly <- midco2$co2_yearly/1000
lowco2$co2_yearly <- lowco2$co2_yearly/1000

legend_ord <- forcats::fct_inorder(c("High","Mid","Low"))

co2plot <- ggplot(highco2[highco2$country == "SE Asia",], aes(x = year, y = co2_cumulative, lty = "High")) +
  geom_line() + 
  geom_line(data = midco2[midco2$country == "SE Asia",], aes(lty = "Mid")) +
  geom_line(data = lowco2[lowco2$country == "SE Asia",], aes(lty = "Low")) +
  scale_linetype_discrete(name = "Scenario", breaks = legend_ord) +
  xlab("Year") + ylab("Cumulative CO2 Emissions since 2016 (Gt)") +
  annotate("label", 2050, highco2$co2_cumulative[highco2$year == 2050 & highco2$country == "SE Asia"],
           vjust = 0.5, label = 'bold("16.61 Gt")', color = "black", size = 3, parse = TRUE) +
  annotate("label", 2050, midco2$co2_cumulative[midco2$year == 2050 & midco2$country == "SE Asia"],
           vjust = 0.5, label = 'bold("12.20 Gt")', color = "black", size = 3, parse = TRUE) +
  annotate("label", 2050, lowco2$co2_cumulative[lowco2$year == 2050 & lowco2$country == "SE Asia"],
           vjust = 0.5, label = 'bold("8.21 Gt")', color = "black", size = 3, parse = TRUE) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1))

pdf("co2plot_seasia.pdf", width = 6, height = 5)
co2plot
dev.off()

midco2$co2_cumulative <- midco2$co2_cumulative*1000
co2plot_2 <- ggplot(highco2_2[highco2_2$country == "SE Asia",], aes(x = year, y = co2_cumulative, lty = "High")) +
  geom_line() + 
  geom_line(data = midco2[midco2$country == "SE Asia",], aes(lty = "Mid")) +
  geom_line(data = lowco2_2[lowco2_2$country == "SE Asia",], aes(lty = "Low")) +
  scale_linetype_discrete(name = "Scenario", breaks = legend_ord) +
  xlab("Year") + ylab("Cumulative CO2 Emissions since 2016 (Gt)") +
  annotate("label", 2050, highco2_2$co2_cumulative[highco2_2$year == 2050 & highco2_2$country == "SE Asia"],
           vjust = 0.5, label = 'bold("19.03 Gt")', color = "black", size = 3, parse = TRUE) +
  annotate("label", 2050, midco2$co2_cumulative[midco2$year == 2050 & midco2$country == "SE Asia"],
           vjust = 0.5, label = 'bold("12.20 Gt")', color = "black", size = 3, parse = TRUE) +
  annotate("label", 2050, lowco2_2$co2_cumulative[lowco2_2$year == 2050 & lowco2_2$country == "SE Asia"],
           vjust = 0.5, label = 'bold("7.16 Gt")', color = "black", size = 3, parse = TRUE) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1))

pdf("co2plot_seasia_2.pdf", width = 6, height = 5)
co2plot_2
dev.off()

co2plot_annual <- ggplot(highco2[highco2$country == "SE Asia",], aes(x = year, y = co2_yearly, lty = "High")) +
  geom_line() + 
  geom_line(data = midco2[midco2$country == "SE Asia",], aes(lty = "Mid")) +
  geom_line(data = lowco2[lowco2$country == "SE Asia",], aes(lty = "Low")) +
  scale_linetype_discrete(name = "Scenario", breaks = legend_ord) +
  xlab("Year") + ylab("Yearly CO2 Emissions (Gt)") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1))
co2plot_annual

pdf("co2plot_annual_seasia.pdf", width = 6, height = 5)
co2plot_annual
dev.off()

######################################## #
#### **** COUNTRY-SPECIFIC PLOTS **** ####
######################################## #

rm(list=ls())

highco2 <- readRDS("highco2.rds")
midco2 <- readRDS("midco2.rds")
lowco2 <- readRDS("lowco2.rds")

highco2$label <- "High"
midco2$label <- "Baseline"
lowco2$label <- "Low"

cdat <- bind_rows(highco2, midco2, lowco2)

co2_seasia <- ggplot(cdat[cdat$country == "SE Asia",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "SE Asia",]) +
  geom_line(data = cdat[cdat$country == "SE Asia",]) +
  xlab("Year") + ylab("SE Asia Total") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_cambodia <- ggplot(cdat[cdat$country == "Cambodia",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Cambodia",]) +
  geom_line(data = cdat[cdat$country == "Cambodia",]) +
  xlab("Year") + ylab("Cambodia") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_indonesia <- ggplot(cdat[cdat$country == "Indonesia",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Indonesia",]) +
  geom_line(data = cdat[cdat$country == "Indonesia",]) +
  xlab("Year") + ylab("Indonesia") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_laos <- ggplot(cdat[cdat$country == "Laos",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Laos",]) +
  geom_line(data = cdat[cdat$country == "Laos",]) +
  xlab("Year") + ylab("Laos") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_malaysia <- ggplot(cdat[cdat$country == "Malaysia",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Malaysia",]) +
  geom_line(data = cdat[cdat$country == "Malaysia",]) +
  xlab("Year") + ylab("Malaysia") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_philippines <- ggplot(cdat[cdat$country == "Philippines",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Philippines",]) +
  geom_line(data = cdat[cdat$country == "Philippines",]) +
  xlab("Year") + ylab("Philippines") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_thailand <- ggplot(cdat[cdat$country == "Thailand",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Thailand",]) +
  geom_line(data = cdat[cdat$country == "Thailand",]) +
  xlab("Year") + ylab("Thailand") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

co2_vietnam <- ggplot(cdat[cdat$country == "Vietnam",], aes(x = year, y = co2_yearly, linetype = label)) +
  geom_line() + 
  geom_line(data = cdat[cdat$country == "Vietnam",]) +
  geom_line(data = cdat[cdat$country == "Vietnam",]) +
  xlab("Year") + ylab("Vietnam") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line.x = element_line(color="black", size = 0.1),
        axis.line.y = element_line(color="black", size = 0.1),
        axis.text = element_text(size = 9)) +
  scale_linetype_discrete(name = "Scenario")

g <- ggarrange(co2_seasia, co2_cambodia, co2_indonesia, co2_laos,
               co2_malaysia, co2_philippines, co2_thailand, co2_vietnam,
               ncol=2, nrow=4, common.legend = TRUE, legend="bottom")

pdf("specificGrid.pdf", width = 6.22, height = 6.36)
g
dev.off()